R Notebook

remotes::install_github("sizespectrum/mizerExperimental")
list<-c("tidyverse", "mizerExperimental", "rfishbase")

lapply(list, require, character.only=T)

## create "not in" operator
'%nin%' = Negate('%in%')

decisions made while exploring:

• use fall biomass
• take a 14 year average of survey observed biomass from start of the time series (the non-fished years)

read in simulated data

survObsBiom <- mskeyrun::simSurveyIndex #atlantisom::read_savedsurvs(d.name, 'survB')
#age_comp_data <- mskeyrun::simSurveyAgeLencomp #atlantisom::read_savedsurvs(d.name, 'survAge') #not using in assessment
len_comp_data <- mskeyrun::simSurveyLencomp #atlantisom::read_savedsurvs(d.name, 'survLen')
#wtage <- atlantisom::read_savedsurvs(d.name, 'survWtage')  #not using in assessment
annage_comp_data <- mskeyrun::simSurveyAgeLencomp #atlantisom::read_savedsurvs(d.name, 'survAnnAge')
annage_wtage <- mskeyrun::simSurveyWtatAge #atlantisom::read_savedsurvs(d.name, 'survAnnWtage')
params <- mskeyrun::simBiolPar
diet <- mskeyrun::simSurveyDietcomp 

#all_diets <- atlantisom::read_savedsurvs(d.name, 'survDiet') #not using in assessment

catchbio_ss <- mskeyrun::simCatchIndex #atlantisom::read_savedfisheries(d.name, 'Catch')
catchlen_ss <- mskeyrun::simFisheryLencomp #atlantisom::read_savedfisheries(d.name, "catchLen")
#fish_age_comp <- #atlantisom::read_savedfisheries(d.name, "catchAge")
fish_annage_comp <- mskeyrun::simFisheryAgecomp #atlantisom::read_savedfisheries(d.name, 'catchAnnAge')
fish_annage_wtage <- mskeyrun::simFisheryWtatAge #atlantisom::read_savedfisheries(d.name, 'catchAnnWtage')

formate for mizer

give species their latin names

sp_params <- data.frame(Name = c("Long_rough_dab",
             "Green_halibut",
             "Mackerel",
             "Haddock",
             "Saithe",
             "Redfish",
             "Blue_whiting",
             "Norwegian_ssh",
             "North_atl_cod",
             "Polar_cod",
             "Capelin"),
    latin = c("Hippoglossoides platessoides",
              "Reinhardtius hippoglossoides",
              "Scomber scombrus",
              "Melanogrammus aeglefinus",
              "Pollachius virens",
              "Sebastes mentella",
              "Micromesistius poutassou",
              "Clupea harengus",
              "Gadus morhua",
              "Boreogadus saida",
              "Mallotus villosus"))

add more species specific info

sp_params <- sp_params %>% left_join(params[,c("Name", "sigma")])

Error in `[.data.frame`(params, , c("Name", "sigma")) : 
  undefined columns selected

get weight at maturity from fishbase

## get mean maturity 
maturity_tbl <- rfishbase::maturity(sp_params$latin)

Joining with `by = join_by(SpecCode)`

## get average for now (can select median, by sex, by location ect)
mean_maturity <- maturity_tbl |>
    group_by(Species) |>
    summarise(age_mat = mean(tm, na.rm = T),
              l_mat = mean(Lm, na.rm = T))

# ## get mean population growth
# popgrowth_tbl <- rfishbase::popgrowth(sp_params$latin)
# popgrowth_tbl <- unique(popgrowth_tbl %>%
#     group_by(Species) %>%
#     summarise(k_vb = median(K)))

sp_params <- sp_params %>% left_join(mean_maturity, by = c("latin" = "Species"))
#sp_params <- sp_params %>% left_join(popgrowth_tbl, by = c("latin" = "Species"))

## calculate w_mat from l_mat
sp_params <-  sp_params %>% 
    mutate(w_mat = a * l_mat ^ b) 

## need age weight relationship for Boreogadus saida who does not have a l_mat
sp_params$w_mat[sp_params$latin == "Boreogadus saida"] <- sp_params$w_max[sp_params$latin == "Boreogadus saida"]*0.25

# ## calculate h 
# sp_params$h <- get_h_default(sp_params)

create generic interaction matrix (use the diet data for this?)

# make a dummy species interaction matrix of full interaction with everything
sp_matrix <- as.data.frame(matrix(rep(1, length(unique(sp_params$species))^2), nrow = (length(unique(sp_params$species))),
                                  dimnames = list(unique(sp_params$species), unique(sp_params$species))))

build mizer model

t <- newMultispeciesParams(species_params = sp_params,
                              #     gear_params = sp_gear,
                                   interaction = sp_matrix, 
                                   initial_effort = 1,
                                   lambda = 2.05, n = 3/4, p = 3/4)

For the species Redfish the value for `w_mat` is not smaller than that of `w_max`. I have corrected that by setting it to about 25% of `w_mat.
Warning: For the following species your value for w_mat and your value for l_mat are not consistent: RedfishWarning: For the following species your value for w_mat and your value for l_mat are not consistent: RedfishWarning: For the following species your value for w_mat and your value for l_mat are not consistent: RedfishNo ks column so calculating from critical feeding level.
Using z0 = z0pre * w_max ^ z0exp for missing z0 values.
Using f0, h, lambda, kappa and the predation kernel to calculate gamma.

run to steady state

ta <- steady(t)

Convergence was achieved in 16.5 years.

calibrate to observed biomass

tb <- calibrateBiomass(ta)
tc <- matchBiomasses(tb)

Warning: For the following species `erepro` has been increased to the smallest possible value: erepro[Green_halibut] = 0.000829; erepro[Saithe] = 0.000417; erepro[Blue_whiting] = 0.00369; erepro[North_atl_cod] = 0.000386; erepro[Polar_cod] = 0.00441; erepro[Capelin] = 0.00497

check plots

plotlySpectra(tc, power = 1, total = T)

plotlySpectra(tc, power = 2, total = T)

plotBiomassVsSpecies(tc)

run to steady state again

td <- tc |>
    calibrateBiomass() |> matchBiomasses() |> matchGrowth() |> steady() |>
    calibrateBiomass() |> matchBiomasses() |> matchGrowth() |> steady() |>
    calibrateBiomass() |> matchBiomasses() |> matchGrowth() |> steady() |>
    calibrateBiomass() |> matchBiomasses() |> matchGrowth() |> steady() |>
    calibrateBiomass() |> matchBiomasses() |> matchGrowth() |> steady() 

Warning: For the following species `erepro` has been increased to the smallest possible value: erepro[Green_halibut] = 0.000829; erepro[Saithe] = 0.000417; erepro[Capelin] = 0.00497Warning: For the following species `erepro` has been increased to the smallest possible value: erepro[Long_rough_dab] = 3.24; erepro[Green_halibut] = 2.39; erepro[Redfish] = 10600; erepro[Blue_whiting] = 0.00551; erepro[Norwegian_ssh] = 0.00387; erepro[Polar_cod] = 0.0873; erepro[Capelin] = 0.055Convergence was achieved in 12 years.
Warning: The following species require an unrealistic reproductive efficiency greater than 1: Long_rough_dab, Green_halibut, RedfishWarning: For the following species `erepro` has been increased to the smallest possible value: erepro[Long_rough_dab] = 20; erepro[Green_halibut] = 13.1; erepro[Redfish] = 61900; erepro[Norwegian_ssh] = 0.012; erepro[North_atl_cod] = 0.000934Convergence was achieved in 16.5 years.
Warning: The following species require an unrealistic reproductive efficiency greater than 1: Long_rough_dab, Green_halibut, RedfishConvergence was achieved in 10.5 years.
Warning: The following species require an unrealistic reproductive efficiency greater than 1: Long_rough_dab, Green_halibut, RedfishConvergence was achieved in 9 years.
Warning: The following species require an unrealistic reproductive efficiency greater than 1: Long_rough_dab, Green_halibut, RedfishConvergence was achieved in 4.5 years.
Warning: The following species require an unrealistic reproductive efficiency greater than 1: Long_rough_dab, Green_halibut, Redfish

summary(td)

An object of class "MizerParams" 
Consumer size spectrum:
    minimum size:   0.001
    maximum size:   25374.6
    no. size bins:  100
Resource size spectrum:
    minimum size:   8.91616e-13
    maximum size:   9.20389
    no. size bins:  175 (221 size bins in total)
Species details:

Fishing gear details:
Gear          Effort  Target species 
 ----------------------------------
knife_edge_gear 1.00   Long_rough_dab, Green_halibut, Mackerel, Haddock, Saithe, Redfish, Blue_whiting, Norwegian_ssh, North_atl_cod, Polar_cod, Capelin 

View(td@species_params)

look at plots

plotBiomassVsSpecies(tb)
plotBiomassVsSpecies(td)

plotlySpectra(tb, power = 1, total = T)
plotlySpectra(tc, power = 1, total = T)
plotlySpectra(td, power = 1, total = T)
plotlySpectra(td, power = 2, total = T)

sim <- steady(td, return_sim = TRUE, t_max = 12, t_per = 0.2)

Convergence was achieved in 0.2 years.

plot(sim)

sim <- steady(td, return_sim = TRUE, t_max = 12, t_per = 0.2)

Convergence was achieved in 0.2 years.
Warning: The following species require an unrealistic reproductive efficiency greater than 1: Long_rough_dab, Green_halibut, Redfish

plot(sim)